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Abstract 

The diagonalization of general mass matrices is a more delicate problem when 
eigenvalue degeneracies exist. In this case, often overlooked in the literature, some 
difficulties arise related to the freedom in the choice of basis in degenerate subspaces. 
Here two simple algorithms are developed to deal with quark and neutrino mass ma- 
trices with arbitrary degeneracies. 

1 Introduction 

In a field theory with fermionic fields ipi, the mass sector of the Lagrangian can be written 
as 

- Cmass = tfiMijlpj + h.c. (1) 

The structure of the mass matrix M is constrained by the symmetries of the Lagrangian: 
charge conservation, baryon and lepton number conservation, etc. We will consider the 
quark and neutrino sectors of the Standard Model with N generations of quarks; ni 
left-handed and ur right-handed neutrinos. The mass terms in the quark sector are 

- C mass = u° Li M? jU % + + h.c, (2) 

where u° L R (d® R ) are the left and right-handed components of the up (down) quark fields, 
and i,j = l,...,N label the different generations of up (down) quarks. The A-dimensional 
mass matrices M u,d are complex and arbitrary. In practice it is convenient to work in a 
quark basis in which M u ' d are diagonal. We make an appropriate change of weak quark 
basis 

ULi = UlijU^j , u Ri = U Rij u° Rj , d Li = U Lij dlj , d Ri = U Rij d^ , (3) 

such that the matrices D u = U"IM U U R , D d = U d M d U R are diagonal, and rewrite the 
mass term as 

- Cmass = ULiDijURj + dliDf^Rj + h.C. (4) 
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Of course this change of basis may affect other terms in the Lagrangian, but we will not 
discuss this here. Thus what we denote here as 'diagonalization' is really a biunitary 
transformation 

U X MU\ = D (5) 

such that D is diagonal (we can also make the matrix elements of D real and positive). 
The matrices U\p, can be calculated []]]] noting that (||) and UiM^U\ = D imply 

UtMM^ul = D 2 , U 2 M ] MU\ = D 2 (6) 

However, this does not determine XJ\ 2 uniquely if D has degenerate diagonal elements, 
and this arbitrariness may cause that D is no longer diagonal: There are matrices U\ t 2 
satisfying (||) that do not necessarily satisfy (||). For instance, let us consider 

«-(?;) p> 

If we try to calculate U 1>2 by using Eqs. (§), we find MM* = M^M = I, so U\ = U 2 = I 
and XJ\MXj\ is not diagonal. This simple example shows us that we need great care to 
deal with the case of degenerate eigenvalues. 
The mass terms in the neutrino sector are (2[ 

- C-mass = 2$LM u lp R + ll.C. (8) 

with 



and My a complex (til + ur) x (til + tir) symmetric matrix. Under a change of basis 

N L = U^ L , N R = U*^ R (10) 
the mass term can be written as 

" Cmass = ^N L D U N R + h.C. (11) 

with D v = UM V U T a diagonal matrix with real positive elements (it will be shown below 
that any complex symmetric matrix M v can be 'diagonalized' with a congruent transfor- 
mation of this type. A proof for the nondegenerate case can be found in Ref. ||). If D is 
nondegenerate we can find the matrix U by diagonalizing M U M^, but if some degeneracies 
are present we can get incorrect results. 
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In this paper we obtain two algorithms for the diagonalization of complex arbitrary 
and complex symmetric matrices, respectively, and we write a package Diagon for their use 
with Mathematica [^]^| This package makes use of Eigensystem to calculate eigenvalues 
and eigenvectors of hermitian and unitary matrices. 



2 Algorithms 



We will first write an algorithm to diagonalize a complex arbitrary matrix by a biunitary 
transformation and then use it to obtain a second one that performs the diagonalization 
of a complex symmetric matrix by a congruent transformation. Note that we could use 
the first algorithm to diagonalize a complex symmetric matrix with two unitary matrices, 
but in general the transformation would not be congruent. 

Let M be a complex arbitrary matrix. It is a well-known fact that we can find two 
unitary matrices Ui, U2 such that U\MU\ = D is a diagonal matrix with real positive 
elements. We will assume for simplicity that D has diagonal elements dj with multiplicity 
nii and ordered from lower to higher values. Let V\, V2 be two matrices that diagonalize 
MM* and M+M respectively, V\MMW^ = D 2 , V 2 M^ MV 2 ] = D 2 . It can be proven that 
Vi = KiUi, where Ki are block-diagonal unitary matrices that commute with D 



( d I, 



D 



"'0 



dil, 



nii 



( id 0) 



Ki 



K 



(i) 



\ dklm k / y 

(the dimension of the blocks is mj). Then, 



ViMVi = K X DK\ = D b 



\ 



(12) 



(13) 



The matrix in the right-hand side of Eq. (13) is a block-diagonal matrix. Each block is 



an unitary matrix times the corresponding eigenvalue of D: Dg 



U) 



more, from DbD^ b 
do = 0, we take Aq 



djK{ j) K ( 2 m . Further- 



D we can compute the eigenvalues di and defining A, = 1/di (if 



1 and D 



(0) 

B 



L m J, 



1/10 



A = diag(A , . . . , X , . . . , X k , . . . , X k ) 



(14) 



1 Although Mathematica has a built-in function SingularValues that performs the diagonalization of 
an arbitrary matrix by a biunitary transformation, it neglects all zero eigenvalues so it is useless for our 
purpose. 
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we find that U 1 =V 1 ,U 2 = AD B V 2 fulfill V x M\j\ = D. 

The diagonalization of a complex symmetric matrix M s by a congruent transformation 
is somewhat more involved. Let U±, U2 be two unitary matrices satisfying UiM s U2 = D, 
with D diagonal with real positive elements (we can calculate U\ 2 as it is shown above). 
As M s is symmetric, U^MgUf = D and these equations imply U2 = KU*, with K a 



block-diagonal unitary matrix that commutes with D (analogous to K{ in Eq. (12)). 
Then, 

UxM s Ul = DK (15) 

We now want to decompose K into two unitary symmetric matrices which commute with 
D. The simplest way to show it is to see that K is a symmetric matrix (all the blocks 
are symmetric except possibly the block corresponding to do = 0, but in this case we 
can redefine U2 without changing U\ and take = I). A complex unitary symmetric 
matrix K can always be written as K = e tS , with S a real symmetric matrix fj, and 



[K, D] = implies [S, D] = 0, so we can rewrite Eq. (15) as 



- l 2\J x M s U{e-^ =D (16) 

However, we need a method to separate the matrix K into two symmetric pieces. First we 
diagonalize K with a real orthogonal matrix O: K = T K^O. It can always be donef] but 
we must check that the subroutine we use to calculate eigenvectors gives real results if it is 
possible (the built-in functions Eigenvectors and Eigensystem do it). We conventionally 
define the square root of K as K 1 / 2 = O t K^ 2 0, with the eigenvectors in the matrix O 
conveniently ordered so that K 1 ^ 2 commutes with D. Finally, U = {K l / 2 )^U\ fulfills 
UMM T = D. 



3 Examples 

Before we present two examples to illustrate the use of the package Diagon, it is impor- 
tant to notice that when the matrix m has degenerate eigenvalues, the set of eigenvectors 
calculated with u=Eigenvectors [m] is not necessarily an unitary matrix, and the prod- 
uct u. Transpose [Conjugate [u]] may have large off-diagonal elements e ~ 10 _1 . This is 
not a roundoff error (typically of order 0(1(T 16 )) as we show in the example below with 
an exact diagonalization. We use DiagonalizeH, a subroutine to diagonalize hermitian 
matrices that uses Eigensystem and the Gram-Schmidt method of orthonormalization in 
each degenerate subspace to yield orthonormal eigenvectors. The following example will 
illustrate this. The numerical results have been obtained running Mathematica 2.2.2 for 
Linux. 

2 Note that S is a real symmetric matrix 
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We write a previously calculated hermitian matrix m with eigenvalues (0,0,1,1) and 
diagonalize it. Note that DiagonalizeH does not return the eigenvectors but their complex 
conjugate. 

In[l] := «Diagon.m 

In [2] := m=l/4 {{1 , (1-1) /Sqrt [2] , (-1+1) /Sqrt [2] , 1} , { (1+1) /Sqrt [2] ,2, 
-1+I,(-1+I)/Sqrt [2]} ,{(-1-1) /Sqrt [2] ,-1-1,2, (-1+1) /Sqrt [2]}, 
{1, (-1-1) /Sqrt [2] , (-1-I)/Sqrt [2] ,3}}; 

In [3] := u=Chop [Eigenvectors [m] ] ; 

In [4] := u.Adj [u] 

Out [4]= {{8, (-2 +2 1) Sqrt [2], 0, 0}, {(-2 - 2 I) Sqrt [2] , 4, 0, 0}, 
8221 221 12 

> {0, 0, -, ( ) Sqrt [2]}, {0, 0, (- + ) Sqrt [2] , — }} 

5 5 5 5 5 5 

In [5] := Simplify [Conjugate [u] . m. Transpose [u] ] 

8 2 2 1 

Out [5]= {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, -, (- + — ) Sqrt [2]}, 

5 5 5 

2 2 1 12 

> {0, 0, ( ) Sqrt [2], — }} 

5 5 5 
In [6] := u2=DiagonalizeH[m] [[2]] ; 
In [7] : = Chop [u2 . Ad j [u2] ] 

Out [7]= {{1., 0, 0, 0}, {0, 1., 0, 0}, {0, 0, 1., 0}, {0, 0, 0, l.» 
In [8] : = Chop [u2 . N [m] . Ad j [u2] ] 

Out [8]= {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 1., 0}, {0, 0, 0, l.» 

The package Diagon (see the Appendix for a complete description) contains the subroutines 
DiagonalizeH, Diagonalize and DiagonalizeS . It has been tested more than 100 times 
on Mathematica 2.2.2 for Linux and Mathematica 2.2 for Windows. We have used the 
instructions listed below with different sets of eigenvalues, obtaining always correct results. 
First we load the package and define the matrices m and ms to test Diagonalize and 
DiagonalizeS respectively. We use 16-digit precision in the input matrices. 

In [9] := SeedRandom[] 

In [10] : = ml=Table [Random [Complex , {-1-1 , 1+I>] , {i , 8} , { j , 8}] ; 

In [11] := m2=Table [Random [Complex, {-1-1,1+1}] ,{i,8},{j ,8}] ; 

In [12] := mal=ml+Adj [ml] ; 

In [13] := ma2=m2+Adj [m2] ; 

In [14] := mul=Chop [Eigenvectors [mal] ] ; 
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In [15] := mu2=Chop [Eigenvectors [ma2] ] ; 

In [16] := m=mul.DiagonalMatrix[{0,0,l,l,l,l,2,2}] .mu2; 

In [17] : = ms=mul . DiagonalMatrix [{0 ,0,1,1,1,1,2, 2>] . Transpose [mul] ; 

In [18] := sol=Diagonalize [m] ; 

In[19] := sol[[l]] 

Out [19]= {0, 0, 1., 1., 1., 1., 2., 2.} 
In [20] := Chop [sol [[2]] .m.Adj [sol [[3]]]] 

Out [20]= {{0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0}, 

> {0, 0, 1., 0, 0, 0, 0, 0>, {0, 0, 0, 1., 0, 0, 0, 0>, 

> {0, 0, 0, 0, 1., 0, 0, 0>, {0, 0, 0, 0, 0, 1., 0, 0>, 

> {0, 0, 0, 0, 0, 0, 2., 0>, {0, 0, 0, 0, 0, 0, 0, 2.}} 
In [21] := sol2=DiagonalizeS [ms] ; 

In [22] := sol2[[l]] 

0ut[22]= {0, 0, 1., 1., 1., 1., 2., 2.} 

In[23] := Chop [sol2 [ [2] ] .ms .Transpose [sol2 [ [2] ]] ] 

Out [23]= {{0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0}, 

> {0, 0, 1., 0, 0, 0, 0, 0>, {0, 0, 0, 1., 0, 0, 0, 0>, 

> {0, 0, 0, 0, 1., 0, 0, 0>, {0, 0, 0, 0, 0, 1., 0, 0>, 

> {0, 0, 0, 0, 0, 0, 2., 0>, {0, 0, 0, 0, 0, 0, 0, 2.}} 
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A Appendix 



(* Title: Diagon *) 

(* Author: J. A. Aguilar-Saavedra *) 
(* Summary: 

This package provides some algorithms for the diagonalization 

of complex arbitrary and complex symmetric matrices 

*) 

(* Version: 1.1 *) 
(* 

This package is available by anonymous ftp at deneb.ugr.es 

in directory pub/packages/ 

*) 

BeginPackage ["Diagon' " , "Linear Algebra' Orthogonalization' "] 

Adj : : usage = 
"Adj [m] gives Transpose [Conjugate [m] ] " 

Ortho : : usage = 

"Ortho [{v,u}] takes a list of eigenvalues v and eigenvectors u and 
returns a list {v,u'} consisting of the eigenvalues v and 
orthonormal eigenvectors u'" 

DiagonalizeH: : usage = 
"DiagonalizeH [m] gives a list {v,u}, with the eigenvalues and 
eigenvectors of the square matrix m, such that 
u.m. Transpose [Conjugate [u] ] =DiagonalMatrix [v] . The list of 
eigenvalues v is ordered by increasing value" 

Diagonalize :: usage = 
"Diagonalize [m] gives a list {v,ul,u2} where v is a row vector, and 
ul,u2 are matrices fulfilling and ul .m. Transpose [Conjugate [u2] ] = 
DiagonalMatrix [v] " 
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DiagonalizeS :: usage = 
" Di agonal izeS [m] gives a list {v,u} where v is a row vector and u is 
a matrix fulfilling u. m. Transpose [u] =DiagonalMatrix [v] " 

Begin [" 'Private' "] 

(* Definitions for internal use only *) 

(* Increase the dimension of a matrix with the identity *) 

IncDiml [m_] : =Transpose [Prepend [Transpose [Prepend [m , Table [0 , 
{Length [m] }] ] ] , Join [{1} , Table [0 , {Length [m] >] ] ] ] 

IncDim [m_ ,nn_] :=Nest [IncDiml ,m,nn] 

(* Pick up a submatrix *) 

DecDim [m_ ,nn_] : =Drop [Transpose [Drop [Transpose [m] ,nn]] ,nn] 

(* Main definitions *) 

Adj [m_] : =Transpose [Conjugate [m] ] ; 

(* Note that the scalar product that works is linear in the 
first variable and antilinear in the second: 
InnerProduct-> (Dot [Conjugate [#2] ,#1]&) *) 

Ortho [sys_] :=Block [{i , j , allvalues , vectors , values , orthvectors} , ( 

allvalues=N [Chop [sys [ [1] ] ] ] ; 

vectors=N [Chop [sys [ [2] ] ] ] ; 

For [i=l , i<=Length [allvalues] , i++ , 

For [j =i+l , j <=Length [allvalues] , j++ , 

If [Chop [allvalues [ [i] ] -allvalues [ [ j ] ] ] ==0 , 

allvalues [ [ j ] ] =allvalues [ [i] ] ] ] ] ; 

values=Union [allvalues] ; 

orthvectors=Flatten [Table [GramSchmidt [vectors [ [Flatten [ 
Position [allvalues , values [ [i] ] ] ] ] ] , 

InnerProduct-> (Dot [Conjugate [#2] , #1] &) ] , {i , Length [values] }],!]; 
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{Sort [allvalues] , Chop [orthvectors] } 

)]; 

DiagonalizeH [m_] : =Conjugate [ 

Drtho[Eigensystem[N[(m+Adj [m])/2] ] ] ] / ; MatrixQ [m] kk 
Length [m] ==Length [Transpose [m] ] 

Diagonal ize [m_] : =Block [{i , soil , sol2 , v,nz,mbox0,mbox,mod,mod2} , ( 
soll=DiagonalizeH[m.Adj [m]] ; 
sol2=DiagonalizeH [Adj [m] .m] ; 
v=Chop [Sqrt [soli [ [1] ] ] ] ; 
nz=Count [v,0] ; 

mboxO=Chop[soll[[2]] .m.Adj [sol2[[2]]]] ; 
mbox=DecDim[mboxO,nz] ; 

mod=Chop [mbox .Transpose [Conjugate [mbox] ] ] ; 

mod2=DiagonalMatrix [Table [1/Abs [Sqrt [mod [ [i , i] ] ] ] , {i , Length [mod] }] ] ; 
{v, Chop [soli [[2]]] , Chop [IncDim [mbox. mod2,nz] . sol2 [ [2] ] ] } 
)] /; MatrixQ [m] kk Length [m] ==Length [Transpose [m] ] kk 
Complement [Flatten [Chop [m] ] ,{0}] !={> 

Diagonalize [m_] : ={Table [0 , {Length [m] }] , IdentityMatrix [Length [m] ] , 
Ident ityMatrix [Length [m] ] >/ ; MatrixQ [m] kk 
Length [m] ==Length [Transpose [m] ] kk 
Complement [Flatten [Chop [m] ] , {0}] =={> 

Diagonal izeS [m_] : =Block [{sol , v,ul ,u2,nz,k0 ,k, sol2 ,rot ,kd} , ( 

sol=Diagonalize [m] ; 

v=Chop[sol[[l]]] ; 

ul=sol[[2]] ; 

u2=sol[[3]] ; 

nz=Count [v,0] ; 

kO=Chop[Conjugate[u2] .Adj [ul]] ; 
k=DecDim[kO,nz] ; 
sol2=0rtho [Eigensystem [k] ] ; 

rot=Sort [Conjugate [sol2[ [2]]] , OrderedQ [{Abs [#2] ,Abs[#l]}]& ] ; 
kd=Chop [rot . k . Adj [rot] ] ; 

{v , Chop [IncDim [Adj [rot] . Sqrt [kd] . rot , nz] . ul] } 

)] /; MatrixQ [m] kk Length [m] ==Length [Transpose [m] ] kk 

Complement [Flatten [Chop [m] ] ,{0}] !={> 
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DiagonalizeS [m_] : ={Table [0 , {Length [m] }] , 
IdentityMatrix [Length [m] ] >/ ; MatrixQ [m] && 
Length [m] ==Length [Transpose [m] ] && 
Complement [Flatten [Chop [m] ] , {0}] =={> 



End[] 



EndPackage [] 
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